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BACKSCATTERING AND NONPARAXIALITY ARREST COLLAPSE OF DAMPED 

NONLINEAR WAVES* 

G. FIBICII 1 *, B. ILAN f S, AND S. TSYNKOV™ 


Abstract. The critical nonlinear Schrodinger equation (NLS) models the propagaiion of intense laser 
light in Kerr media. This equation is derived from the more comprehensive nonlinear Helmholtz equation 
(NLH) by employing the paraxial approximation and neglecting the backseat tered waves. It is known that 
if the input power of the laser beam (i.e., L 2 norm of the initial solution) is sufficiently high, then the 
NLS model predicts that the beam will self-focus to a point (i.e., collapse) at a finite propagation distance. 
Mathematically, this behavior corresponds to the formation of a singularity in the solution of the NLS. A key 
question which has been open for many years is whether the solution to the NLH, i.e., the ‘‘parent' equation, 
may nonetheless exist and remain regular everywhere, in particular for those initial conditions (input powers) 
that lead to blowup in the NLS. In the current study wo address this question by introducing linear damping 
into both models and subsequently comparing the numerical solutions of the damped NLH (boundary-value 
problem) with the corresponding solutions of the* damped NLS (initial-value problem). Linear damping 
is introduced in much the same way as done when analyzing the classical constant-coefficient Helmholtz 
equation using the limiting absorption principle. Numerically, we have found that it provides a very efficient 
tool for controlling the solutions of both the NLH and NLS. In particular, we have been able to identify 
initial conditions for which the NLS solution does become singular, whereas the NLH solution still remains 
regular everywhere. We believe that our finding of a larger domain of existence for the NLH than that for 
the NLS is accounted for by precisely those mechanisms that have been neglected when deriving the NLS 
from the NLH, i.e., nonparaxiality and backscattering. 

Key words. Kerr medium, nonlinear wave propagation, self-focusing, singularity formation, linear 
damping, limiting absorption, two-way ABCs 

Subject classification. Applied and Numerical Mathematics 
1. Introduction. The focusing critical nonlinear Schrodinger equation (NLS) 

i^z(z, x) + A_l 0 + |0| 4/rf 0 = 0, 0(0, x) == 0o(x), (LI) 

where x £ R d and Aj_ = d X]X} H h d XdXd , arises in a variety of physical contexts. Of foremost interest is 

the ease d = 2, which corresponds to the propagation of intense laser beams in Kerr media. In this ease, z is 
the axial coordinate in the direction of propagation, x = (x.y) are the spatial coordinates in the transverse 
plane. Ax = d xx + d yy is the diffraction term (transverse Laplacian). and |0| 2 0 describes the nonlinear 
polarization of the Kerr medium. It is well known that solutions to t.he critical NLS (1.1) can self-focus 
and eventually collapse, i.e., become singular, at a finite propagation distance, provided that their initial 
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power A r (0) — f |^o| 2 dx exceeds a threshold power Ay, whose value depends only on the dimension d [7,28]. 
Since, however, physical quantities do not become infinite, and since in experiments laser beams continue to 
propagate beyond the NLS blowup point, the question arises as to what specific physical mechanism(s) among 
those that have been neglected when deriving the NLS from the Maxwell's equations, actually arrest(s) the 
collapse. We recall that the final stage in the derivation of the NLS is to disregard the backscattering and 
apply the paraxial approximation (see Section 2.2) to the critical nonlinear Helmholtz equation (NLH) 

±E(z, x) + k% (l + e|£| 4/d ) E = 0, A ee d zz + A ± . (1.2) 

where k 0 is the linear wavenumber and the extent of nonlinearity is measured by the quantity c = 4co^ 2 , 
where t? 2 is the Kerr coefficient, see, e.g., [3, 19]. Therefore, it is natural to ask whether going back from 
the NLS to the NLH, i.e., adding nonparaxiality and backscattering, is sufficient to guarantee existence of 
the solution with no singularities. In other words, for a given initial condition that leads to blowup in the 
critical NLS. does the NLH (always) have a solution that remains regular everywhere? 

The foregoing question has been open for many years. In his celebrated 19G5 paper in Physical Review 
Letters [15], which was the first publication in the literature to predict that the solutions to the critical 
NLS can become singular, Kelley was careful to note that the paraxial approximation, and hence the 
entire NLS model, breaks down near the singularity. Frit and Fleck [4] were the first to demonstrate 
that nonparaxiality of the beam can arrest the blowup, by showing numerically that the initial conditions 
that lead to singularity formation in the NLS, result in focusing-defocusing oscillations in the NLH. In these 
simulations, however, they did not solve a true boundary- value problem for the NLH. Instead, they solved an 
initial- value problem for a “modified” NLH that only describes the right-propagating wave (while introducing 
several additional assumptions along the way). Akhmediev and collaborators [1,2] analyzed an initial- value 
problem for a different “modified” NLH; their numerical simulations also suggested that nonparaxiality 
arrests the singularity formation. Both numerical approaches [4] and [1,2], however, did not account for the 
effect of backscattering. Fibich [5] applied asymptotic analysis to derive an ODE in z for self-focusing in 
the presence of small nonparaxiality. His analysis suggests that nonparaxiality indeed arrests the singularity 
formation, resulting instead in decaying focusing-defocusing oscillations. However, backscattering effects 
were neglected in this asymptotic analysis. 

The aforementioned studies [1,2,4, 5, 15] have prompted a general belief that nonparaxiality arrests the 
collapse. However, no rigorous proof of global existence for the NLH has ever been provided. Moreover, all the 
simulations in the above studies neglected the backscattering and considered only the forward-propagating 
field. The first numerical solutions of the NLH as a true boundary value problem, with backscattering 
effects fully included, have been obtained by Fibich and Tsynkov in [12], using a high-order discretization 
supplemented by a new two-way artificial boundary condition (ABC). The simulations in [12] were performed 
for the values of the input power of up to 90% of the threshold Ay, and they have captured the mild self- 
focusing of the corresponding solutions. In a subsequent paper [10], we have corroborated experimentally 
the prediction of the asymptotic analysis that the magnitude of the backscattered signal scales quadratically 
with the nonparaxiality parameter / (see Section 2.2), and that the computed NLH solutions converge to 
the corresponding NLS solutions as / goes to zero. 

The numerical methodology of [12] was obviously not free of limitations of its own. Foremost, we could 
not obtain converging solutions for initial powers equal to or higher than the critical value Ay. In [12], we 
have considered initial powers of only up to 90%. of N c \ in the current paper we computed the NLH solutions 
for up to A' (0) = 0.99 A' c (see Section 4). In the course of these simulations we have noticed that as A r (0) 


2 



approaches the critical power from below, the convergence rate of the iterations slows down noticeably. This 
makes the simulations for higher subcritieal values of A'(0) (0.99AV < A r {0) < Ay) difficult to conduct, 
although it is reasonable to assume that the NLH solution will converge 1 for input powers all the way up to 
N c . However, for the input power A r (0) exactly equal to Ay the converge nee of nonlinear iterations of [12] is 
lost (see Section 4). 

The aforementioned slowdown of convergence for input powers slightly below N c should be attributed 
to either deficiencies of the method, or to insufficient computational resources, or to both. As concerns 
the iteration method of [12] itself, it is the most straightforward approach based on simply freezing the 
nonlinearity; most likely, it can be improved or replaced by a more advanced technique, and we plan on 
looking into this issue in the future. As for the computer resources requirements, they are determined by the 
size of the computational domain, which should be sufficiently large so that to meet the condition of near- 
linear propagation in the far field, see [12]; and by the grid size, which should be sufficiently fine to resolve a 
given wave length and the sharp near-blowup profile. These requirements become more stringent for higher 
input powers, which decay at larger distances and/or undergo stronger focusing. In other tvords, the higher 
the input power the larger the domain and/or the finer the grid that one reeds to use in order to maintain the 
same solution quality and/or convergence rate. In our previous simulations w'e have, indeed, seen examples 
of diverging NLH solutions with subcritieal input powders which converged on a larger computational domain 
and/or at a finer resolution. It is still unclear, however, whether having more computer resources and/or a 
better nonlinear iteration scheme will allow one to solve the NLH for initial conditions that lead to collapse 
in the NLS, or whether the convergence breakdown at A r (0) > A r r is an indication of the loss of solvability 
of the NLH, or loss of regularity of the solution. 

As such, in the current paper we explore an alternative approach t o the issue of solving the NLH in 
the blow up regime of the NLS, by considering the linearly damped NLH and the corresponding linearly 
damped NLS. The addition of linear damping is not an ad hoc procedure. Indeed, an electromagnetic wave 
is always partially absorbed by the medium through which it propagates, an effect neglected in either the 
original undamped NLH or NLS, both of wdiich model the propagation under “ideal transparency.” A 
mathematical motivation to add linear damping comes from the so-called limiting absorption principle that 
is used for identifying the unique solutions of the linear Helmholtz equation, see, e.g., [27]. It is known that 
the classical constant-coefficient homogeneous Helmholtz equation 

+ = 0 (1.3a) 

has non-trivial solutions on the entire space even in the class of functions that vanish at infinity, which 
obviously amounts to non-uniqueness. To fix the problem, the additional Sommerfeld boundary conditions 
need to be introduced at infinity that basically distinguish between the incoming and outgoing waves. On 
the other hand, wdien a complex absorption coefficient is added, the new damped equation 

AE + *g(l + iS)E = 0 (1.3b) 

has only trivial solution. Consequently, its inhomogeneous counterpar will be uniquely solvable for any 
compactly-supported right-hand side in rather wide classes of functions, such as tempered distributions, 
see [27]. Moreover, when S — > ±0, the unique solution of the inhomogeneous damped equation will converge 
uniformly on the entire space to the solution of the respective undamped equation that corresponds to either 
the radiation of waves toward infinity (outgoing waves), or conversely, the incidence of weaves from infinity 
(incoming waves), where the distinction is rendered by the sign of 6 . This, in particular, implies that if 
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wo decide to keep a small but finite damping in the equation, we may expect its solution to be uniformly 
close to the solution of the undamped equation that is driven by the same source terms and is composed of 
either only outgoing or only incoming waves in the far field. The latter consideration is especially important 
in the context of our iteration algorithm, see Section 3 and [12] for detail, which basically reduces to a 
repeated solution of the constant-coefficient. Helmholtz equation driven by a variety of compactly supported 
right-hand sides and subject to the radiation boundary conditions in the far field. 

Solving the damped NLH numerically as a true boundary value problem required only minor changes 
in the algorithm of [12] for the undamped NLH, which are described in Section 3. At the same time, the 
addition of damping allows us to better control the solution. In particular, damping decreases the solution 
magnitude in the far field, which is a key requirement for the validity of the artificial boundary conditions 
(ABCs) of [12]. As a result, we have been able, to consider initial conditions with the powers well above N c . 

Let us recall that for a given initial condition that leads to the blowup in the undamped critical NLS. 
there is a threshold value <S t s h of the damping parameter 8 such that if 8 > then linear damping arrests 
the collapse, whereas when 8 < 8f h the solution of the NLS blows up, see [G]. 1 In the numerical simulations 
of the damped NLH reported hereafter we found a similar threshold value such that for 6 > ^ the 
solution exists and is regular everywhere, whereas when 8 < <5^ the iteration scheme diverges. As has been 
mentioned, in the latter case it is not clear whether the divergence indicates that there is no solution to 
the NLH, or that our computational resources art 1 insufficient (or the iteration scheme is suboptimal) to 
calculate the solution. Therefore, we can conclude that the actual (analytical) threshold value Jjj,, such that 
regular solutions to the NLH exist for all 8 > Jjj lT is less or equal than the computed threshold which is 
determined from the simulations, i.e., that 0 < fill < fill 

The main result of the current study is that 

r H ^ rS 

”lh < ^th- 
in other words, for a given initial condition that leads to the blowup in the undamped NLS, there is an entire 
range of values for the damping coefficient: < 8 < £ t s h , for which the damped NLS solution will blow 

up, but the NLH solution will be regular everywhere*. Therefore, we can conclude that nonparaxiality and 
backscattering arrest the collapse when the damping parameter is in the range 8^ h < 8 < 8f h . Whether NLH 
solutions exist for infinitely small linear damping as well, i.e., in the limit 8 — > 0, is the question that yet 
remains to be answered. We believe, however, that this question should be considerably easier to address, 
both numerically and analytically, than the question of solvability of the original undamped NLH. 

2. Formulation of the Problem. 

2.1. The Nonlinear Helmholtz Equation. A typical setup for the propagation of electromagnetic 
waves in a Kerr medium is shown in Figure 2.1. An incoming laser beam with known characteristics impinges 
normally on the planar interface z = 0 between the linear and the nonlinear medium. The electric field E — 
E(z,x) is governed by the nonlinear Helmholtz equation (1.2). For simplicity, we consider the cylindrically- 
svmmetric case 2 where E — E(z,r) and r = + . . . + x 2 (i . The nonlinear medium occupies the semi-space 

z > 0 (see Figure 2.1). Consequently, the NLH (1.2) has to be supplemented by boundary conditions at z = 0 

1 Self-focusing in the critical NLS is highly sensitive to the effect of small perturbations. Some perturbations can arrest the 
collapse even if they are initially infinitesimally small [11]. In contrast, an infinitesimally small linear damping does not arrest 
the collapse, and its sufficient amount must be present to regularize the solution. 

2 This assumption is quite reasonable, since even when the initial conditions of the NLS are not cylindricallv-symmetric, 
near the singularity the solution becomes cylindrically-symmetric [8]. 
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KlG. 2.1. Schematic of propagation of waves in Kerr media. 


arid z — > -foe. We require that as 2 — > + 00 , E have no left-traveling components and that the propagation 
be diffraction-dominated with the field amplitude decaying to zero, i.e., lim max \E(z,r)\ = 0, which also 

; -4x 0<r<x 

means that the nonlinear wavenumber k 2 = k.Q (l + f \E\ A ^ d ) approaches its linear limit: lim k 2 = In 

other words, at large z' s the solution should be a linear superposit ion of right-traveling waves. Since the 
actual numerical simulation is carried out on a truncated domain 0 < 2 < 2 max (Figure 2.1). the desired 
behavior of the solution as z — > -hoc has to be captured by a far- field artificial boundary condition (ABC) 
at the artificial boundary z = 2 ma x- This boundary condition should guarantee a reflectionless propagation 
of all the waves traveling towards z = -hoc. Often, boundary conditions designed to ensure the transparency 
of the outer boundary to the outgoing waves are called radiation boundary conditions [24]. 

The situation is more complex at the interface z = 0. where the 1 total field £(0, r) is composed of a given 
incoming (right-traveling) component E\ nc (Q,r) and an unknown backscattered (left-traveling) component 

£scat,(0. r), i.e., 


£(0, r) = £i nc (0, r) -h J5 sca t(0,r). 

As such, the boundary condition at . z = 0 has to guarantee the reflect ionless propagation of any left- 
traveling wave through the interface and at the same time be able to correctly prescribe the incoming signal. 
Implementation of such a two-way ADC was first carried out in [12] for the undamped NLS, and is extended 
to the damped case in Section 3.3. 

Finally, the electric field vanishes as r — > -hex:. In practice, we truncate the domain at some large but 
finite r, nax and require that £(z,r max ) = 0. 

2.2. Paraxial Approximation and the Nonlinear Schrodingei Equation. We first introduce the 
dimensionless quantities f. 2 , and t/’ as 

f=~.. 2 = E = ^ s {f.r 2 X)- <l t\>{s,r), ( 2 . 1 ) 

r 0 
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where ro is the transverse width of the input beam and Lpy ~ k 0 r % is the diffraction length. Then, by 
substituting the quantities (2.1) into the NLH (1.2) and dropping the tildes, we obtain 

in>: + A±V' + W 4/< ty = -4 f 2 U>zz, (2.2) 


where / = 1 /r^ko = X/2nro is the nonparaxiality parameter. 

The standard derivation of the NLS is motivated by the observation that / « 1, since typically A <$: r 0 . 
This suggests that one can neglect the y> zz term, i.e., apply the paraxial approximation , and obtain the 
nonlinear Schrodinger equation 

iip z (z,r) + = 0, (2.3) 

which is the same as the previously introduced equation (1.1), except that in (2.3) we use r instead of x for 
simplicity. The NLS (2.3) is supplemented by the initial condition at z — 0 : 

0(0, r) = E- mc {Q,r). 


Subsequently, it needs to be integrated by a “time” -marching algorithm, where the direction of propagation 2 
plays the role of time. We reemphasize that hackscattering effects are not taken into account by the NLS (2.3). 
Indeed, once (2.3) is solved, the overall solution, according to (2.1), is the slowly varying amplitude t times 
the forward propagating oscillatory component e tk()Z . 

2.3. Linear Damping. When damping, i.e., linear absorption, is included, the NLH (1.2) becomes 

A E(z, x) + kU 1 + iS + f|£| 4/d ) E = 0, (2.4) 

where ko is the (real part of the) wavenumber, 

O(rag) 

and no is the linear index of refraction of the medium. The corresponding NLS (2.3) becomes (see eq. (2.1)) 

ixpz + 4- \y>\ 4 / d ip -I- irokodip = 0. (2.5) 


By definition, optical transparency of the medium means that the damping is small. For example, for water 
in the visible regime [14], 


3(n.p) 


10 “ 7 . 


Having small physical values of damping also agrees well with the mathematical reasoning behind the limiting 
absorption principle. As indicated in Section 1 (see, e.g., [27] for detail), for a classical constant-coefficient 
Helmholtz operator of (1.3a), the introduction of a small complex absorption coefficient of the appropriate 
sign [as in (1.3b)] implies that there will be a unique solution for any compactly supported excitation, and 
that this solution will be uniformly close in the entire space ! d+l to the solution of the corresponding 
undamped linear Helmholtz equation driven by the same sources and subject to the radiation boundary 
conditions in the far field. In the following Section 3, we show that for the formulation analyzed in this 
paper the proper sign of 6 is positive. 

As we have noted before, the physical case that corresponds to the propagation of laser beams in bulk 
Kerr media is d = 2. However, in order to reduce the complexity of the computations we rather consider a 
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simpler case d — 1, as was previously done in [12]. Thus, the damped XLH for E = E(z,r) and the damped 
NLS for v? = ty(z,r) that are solved numerically in this study are 

E zz (z,r) + Err + A'q ( 1 + is + f|£| 4 )£ = 0, (2.6) 


and 


iii> z (z, r) + tVr + ir^k^Sv + |V>| 4 V = 0. 


(2.7) 


respectively. 

3. Numerical Methods. The damped XLH (2.G) is solved using fourth-order finite differences. The 
methodology of solution is outlined below in Section 3.1; it is similar to the one that we have introduced in 
our previous work [12] for solving the undamped NLH. The choice of a higher-order method is motivated 
primarily by the necessity to resolve a small-scale phenomenon of backscattering at the background of the 
forward propagating waves. The damped NLS (2.7) is also solved by a fourth-order scheme; it is natural to 
expect that this will leave less room for potential purely numerical discrepancies between the two techniques 
and as such, will allow for a more accurate comparison. Besides, it is generally known that higher-order 
methods provide for a better resolution of waves. 

3.1. Discretization of the NLH and Solution Methodology. We use a conventional fourth-order 
central-difference discretization for the Laplacian A = d zz + d rr of (2.6) in so doing the stencil is five-node 
wide in each coordinate direction. As the equation is nonlinear, we implement a nested iteration scheme. 
On the outer loop, we freeze the nonlinearity, i.e., consider the coefficient k 2 = &q (l -f iS T- f|£ , | 4 ) as a 
given function of the coordinates 2 and r, which is actually obtained by taking the quantity \E\ 4 from the 
previous iteration, see (2.6). This way we arrive at a linear equation with variable coefficients. The' latter is 
also solved by iterations on the inner loop of the nested scheme. Here, we leave the entire varying part of 
the equation, which is proportional to e, on the lower level, and on the 1 upper level need to invert only the 
constant- coefficient linear damped Helmholtz operator A + k^(l + iS)I (cf. equation (1.3b)). 

Formally,' our iteration scheme resembles the fixed-point approach, however, no rigorous convergence 
theory is available yet. and the convergence is assessed experimentally. The advantage of using these' nested 
iterations is that first, the method eventually reduces to the repeated solution of one and the same linear 
constant coefficient equation driven by different source terms, which can be done efficiently at the discrete 
level. Second, the radiation boundary conditions at 2 = z max and the two-way ABCs at z = 0 . see Figure 2.1, 
are most convenient to set on the upper time level of the iteration scheme already for the linear constant- 
coefficient operator. 

To solve the linear constant-coefficient damped discrete Helmholtz equation 

A ih) E + k%(l + i6)E = g, (3.1) 


where g is the right-hand side generated on the previous iteration, we first separate the variables by 
implementing the discrete Fourier transform in the transverse direction r; the boundary conditions are 
symmetry at r = 0 and zero Dirichlet at r = r max (see Section 2.1). This yields a collection of fourth-order 
one-dimensional finite- difference equations (grid index n corresponds to the continuous variable z): 


E n _ 2 + 16£„_i - 30 E n 4- 1GE„ +1 - E n+2 
\2h 2 


+ (At 0 (l + IS) A m )E n — g n 


(3.2) 



parameterized by the dual Fourier variable A m , the latter is defined by formula (29) of [12]. Each equation 
(3.2) needs to be solved independently. 3 The two-way and radiation ABCs at z = 0 and ~ — 2 niax , 
respectively, for the discrete equation (3.1) are set in the Fourier space, i.e., individually for each one- 
dimensional equation (3.2). This is done by first identifying the linearly-independent eigen- modes for the 
homogeneous version of this equation. It is important, to note that even though the original differential 
equation is of the second order, we are using its fourth-order approximation and as such, each homogeneous 
discrete one-dirnensional equation of type (3.2) has four' linearly independent solutions. These solutions are 
q”, n , </.?, and q 2 r \ see [12], where q x , 1 /q\, q 2 , and 1 /q 2 are roots of the characteristic algebraic equation 

-1 + I6q + (I2h 2 z (kl(l + iS) - A m ) - 30)g 2 + 1 6q 2 - q 4 = 0. (3.3) 

3.2. Roots of the Characteristic Equation. It is indeed easy to see that equation (3.3) has two 
pairs of mutually inverse roots. We first notice that this equation originates from a central-difference, i.e., 
symmetric, discretization (3.2). As such, if q is a root, then q~ x is obviously a root as well, which can be 
verified by direct substitution. Then, to actually find the roots we rewrite the polynomial on the left-hand 
side of (3.3) as 

(9 - 9i)(<7 - - Q2)(q - (l2 X ) 

= - 1 + (di + d 2 )q - (2 + d { d 2 )q 2 + ( d\ + d 2 )q 3 - q 4 , 


where 


di = <7i + <zr\ d 2 = q -2 + 9 2 -\ 

and match the coefficients. In so doing, we obtain 

d x + d 2 = 16, -2 - d x d 2 - I2hi(kl(l + iS) - A ra ) - 30, (3.4) 

so that each pair of roots: q \ , q f 1 and q 2 , q 2 } , can be found by solving the corresponding quadratic equation: 

q 2 — d x q -1-1 = 0 (3.5a) 

or 

q 2 — d 2 q +1=0, (3.5b) 

while the coefficients d\ and d 2 are, in turn, determined by solving quadratic equations (3.4). 

At this stage, the key difference between the current analysis for the damped equation and the previous 
analysis for the undamped equation of [12] needs to be emphasized. As shown in [12], when 5 = 0 the first 
pair of solutions of the homogeneous equation (3.2), q [* arid gf”, approximates the genuine “longitudinal,” 
i.e., ^-aligned, modes of the undamped homogeneous differential equation (1.3a): 

E x = e ikeZ , and E 2 = e’ ik * z , (3.6) 

respectively. The functions E\ - E\{z) and E 2 — (z) are two linearlv-independent solutions of the ODE 

E zz + (hi -X)E = 0 (3.7) 

3 Note, the discrete equations (3.1) and (3.2) are very similar to the corresponding discrete equations studied in [12] except 
that previously we had no damping. 
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obtained by Fourier transforming equation (1.3a) with respect to r; \ is the dual variable. In formulae 
(3.6), we have denoted k c = y/kfi - A, and a particular branch of the square root that we always take is 
\/ pe i0 = p l /' 2 e i0 / 2 . The two continuous modes (3.6) may be either traveling or evanescent waves depending 
on whether the real quantity k 2 - (k$ - A) is positive or negative, or in other words, whether the dual Fourier 
variable A is less or greater than fcg. To demonstrate' the aforementioned approximation property for the 
undamped (6 = 0) discretization (3.2), we re-define k r = y/ti 2 - A„~, introduce a = fizk and show in [12] 
that if a > 0 then q] and qf 1 are complex conjugate roots of the characteristic equation (3.3). Both these 
roots have unit magnitude: |r/i| = \qf l \ = 1, which indicates that q" and q x T1 are pure discrete traveling 
waves. Moreover, if a 1 then (see [12]) 

<71 = ^ v ' ,£ + 0({k c h z )% f/f 1 = f + 0((k c l h f). (3.8) 

Equalities (3.8) imply that in the undamped case S = 0, q” is a discrete counterpart of the right-traveling 
wave E[, and q^ n is a discrete counterpart of the left -traveling wave E ? ; the approximation is obviously 
fourth-order accurate because on the grid z n = h z ii. If « < 0 and still S == 0, then we again show in [12] that 
|f/i | < 1 and \qf ] \ > 1, which indicates that q[ l is a right-evanescent wav-’ and qf n is a left-evanescent wave. 

The situation changes drastically with the introduction of damping. In contradistinction to the 
undamped case, when 6 ^ 0 the homogeneous differential equation no longer has pure propagating, i.e., 
constant-amplitude, longitudinal modes. Indeed, by Fourier transforming equation (1.3b) in the r direction, 
we arrive at the family of ODEs 

^« + (*g(l+W)-A)£? = 0 (3.9) 


parameterized by the dual variable A. Each of the equations (3.9) has two linearly independent solutions: 

(3.10) 


£, = e iW*?+U'Z» = e li 


■4 


+ 1 - 4 * <5 


t — p- i: y/k?+ ik Q d — ( 


ik r z t/ 1 H- # 


Clearly, the second equality in each formula (3.10) is valid only if k c ^ 0 Formulae (3.10) show that as long 
as 8 ^ 0 there will always be a nontrivial real part in each exponent. Consequently, the amplitudes of the 
waves (3.10) will always decrease or increase exponentially for z — > ix. In particular, if we analyze the 
traveling waves regime of the undamped equation, i.e., the case of small A: k$ — A > 0, and additionally 
assume that \S\ 1, then formulae (3.10) yield (cf. formulae (3.6)): 


A{ damped) 

^ ( , ikr ~ 

_ f ik rZ-i^6z _ ^undamped) _ Sz 

2 ~r( damped) ^ 

-ikrzf 1 + ii 

t e v ‘ k r , 

) _ e -ik r z+^Sz _ ^(un, lamped) 


(3.11) 


Since we identify £‘ undBmped) = e ik ' z of (3.6) as the right-traveling wave, and £< undamppd) = e~ ik ' : of (3.6) 
as the left traveling wave, we can conclude that to have the propagation toward infinity (i.e., the radiation 
of waves) accompanied by the decay of the amplitude (as opposed to growth with no bound), we have to 
take positive values of the damping factor: S > 0 (cf. Section 1). In this case, the amplitude of £< (lam P pd) 
will decay exponentially for z — > +oc (propagation to the right), and the amplitude of £^ dam P pd) will decay 
exponentially for 2 — > -oc (propagation to the left). As one can easily sec' from (3.11), the rate of decay is 
controlled by the value of 6 . 
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In connection to the aforementioned exponential behavior of the longitudinal modes, a more general 
fact is also worth mentioning. The full Fourier symbol of the undamped operator of (1.3a) obviously has 
real roots on the dual plane, these roots occupy the entire circle of radius k 0 centered at the origin. In 
contradistinction to that, the symbol of the damped operator of (1.3b) does not have real roots on the dual 
plane. As shown in [20], the damped operator will therefore have an exponentially decaying fundamental 
solution. In practical terms it means that the outgoing waves governed by the damped Helmholtz equation 
will decay exponentially toward infinity in all directions. For comparison we remind that the fundamental 
solution of the undamped operator is given by a zero order Hankel function, which only decays at infinity 
as the inverse square root of the distance from the origin. 

To establish the properties of the propagating modes for the discretization (3.2) in the presence of 
damping, and to demonstrate similarities to the continuous damped case, we first introduce and prove 
Proposition 3.1. The characteristic equation (3.3) for <5/0 does not have roots with unit magnitude. 
Proof. Let us assume the opposite: There exists a unit magnitude root q = e tB to the algebraic 
characteristic equation (3.3). Then, 

- 1 + 16c" + (12/i? (A:q ( 1 + iS) - A m ) ~ 30)c 2 ^ + ICC 3 ** - e 4 ^ 

= [-e- 2 " + 16e"** + (12h*(Jfcg( 1 + iS) - X m ) - 30) + 16e 10 - e 2ie ] • e 2 ‘* 

= [—2 cos(20) + 32 cos 0+ (12h 2 (^(l + *<$) - X m ) -30)] • e 2 ^ = 0. 

As e 2ie / 0, the expression in rectangular brackets has to be equal to zero. Since the only imaginary 
contribution to this expression is proportional to <5, we conclude that it is only possible when 8 — 0. □ 

Proposition 3.1 implies that similarly to the continuous case, there will be no constant-amplitude 
solutions to the homogeneous counterpart of the discrete equation (3.2). Each of the four corresponding 
modes: q\\ gf n , q 2 , and q 2 n , will exponentially decrease in one direction and exponentially increase in 
the opposite direction. In particular, if we assume as before that a <§C 1 in the undamped traveling waves 
regime, 4 and in addition let <5 <C 1, then solving first equations (3.4) for d\ , then equation (3.5a) for <71 and 
q x 1 , and finally using the Taylor expansion, we obtain (cf. formula (3.8)) 


<71 =e ik ' tl >-^ Sh ‘ +0 k c h z 


q -\ + q( kchz 


Equalities (3.12) mean that the damped discrete traveling waves q™ and q(~ n approximate the damped 
continuous waves (3.11) with the fourth order of accuracy. This result is obviously similar to the one 
obtained in the undamped case, see formulae (3.8). 

As of yet, our discussion has focused on the first pair of roots q\ and gf 1 of the characteristic equation 
(3.3), because these roots correspond to the genuine modes of the original differential equation. The second 
pair of roots q 2 and q% 1 is obtained by solving equations (3.4) for do and subsequently solving equation 
(3.5b). The corresponding pair of solutions q " and q^ 71 is, of course, a pure numerical artifact. In [12] we 
have' shown that for <5 = 0 the roots q\ and q(~ l cannot have unit magnitude: |g 2 | < 1 and | q 2 l \ > P which 
means that the waves q” and q 2 n are always evanescent. In the damped case, Proposition 3.1 implies that 



is small and k c ~ ko. 


10 


these waves will remain evanescent as well. The presence of the second pair of waves, however, implies that 
the discrete equation requires two more boundary conditions compared 10 the original differential equation. 

In Section 1 , we have outlined a general two-fold motivation behind the introduction of damping into 
the Helmholtz equation. One part was coming from physics because absorption by the medium always 
accompanies the propagation of electromagnetic waves in real-life settings. Moreover, from the standpoint 
of mathematics the introduction of damping helps select a unique solution using the limiting absorption 
principle. Besides these two key reasons, the presence of damping in the equation also affects positively the 
properties of the numerical algorithm. 

First, having no roots of unit magnitude presents a significant advantage from the viewpoint of numerical 
stability . In this case, every discrete system (3.2) supplemented by the boundary conditions that are discussed 
below in Section 3.3, will be well-posed in the classical sense of [13.21]. In contradistinction to that, in the 
original undamped case existence of the roots with unit magnitude may, generally speaking, cause a weak 
polynomial growth of the error when the grid size is refined, although no major exponential instability will 
be possible. 

Second, we remind that the original formulation of the problem requi es that E(z , r) vanish as \r\ — » oc. 
Instead, when solving the problem numerically we set E(z , r) = 0 at a laige but still finite distance r = r max . 
Of course, we expect that on some fixed bounded region of interest located next to the axis of the propagating 
beam our solution will converge to the original infinite-domain solution with the increase of r max . A general 
methodology for solving infinite-domain problems based on a similar idea was first introduced and studied 
in [22,23,25,26] in the context of fluid flow. It was shown, in particular, t .iat one may obtain the convergence 
rate inversely proportional to the square of the domain size (i.e., ~ l/ry ax using our particular notations). 
Besides, for a specific example that involves the Laplace equation that transforms into a Yukawa equation 
by introducing small “dissipation,” Mishkov and Ryaben’kii have shown in [18] that one may expect a much 
faster convergence of the damped solution to the undamped one on a fixed-size domain rather than on the 
original unbounded domain. Even though the formulation of the problem in [18] is not quite the same as the 
one analyzed here, there are still similarities that allow us to consider the results of [18] as another argument 
toward using the damped equation. 

3.3. Boundary Conditions. Apart from the foregoing key difference in the properties of the roots of 
equation (3.3) in the undamped and damped case, see Section 3.2. the algorithm for solving the damped NLH 
remains basically the same as the undamped algorithm of [12]. Each equation (3.2) needs to be supplemented 
by the radiation boundary conditions at 2 = z mSLX and two-way ABC's at z — 0. 

The radiation boundary conditions are constructed by requiring that on the right boundary z — z m ax the 
solution of (3.2) be composed of only the waves that propagate/ decay to the right, i.e., E n = c\q\ l + c^q" • 
The selection is rendered by the so-called one-way discrete Helmholtz equation [12], which is a a linear 
homogeneous relation that defines the span of all the appropriate modes. Specifically, let us consider equation 
( 3 . 2 ) on the grid n = 0, 1, . . . , N — 1, Ay and assume that the right-hand side g n is small and can therefore 
be neglected near the right boundary n = Ay i.e., that the propagation is almost linear in the far field. 
Then, we require that the vector [£yv_ 3 , En- 2 , En~\ , En] t be a linear combination of the two vectors: 
[q\ ~'^q\'~ 2 ,q?~ l ,q]'] r and [q£ “ 3 , <z ^ -2 , ^ “ 1 , q^] 1 > which obviously translates into 

En- 3 £/v -2 En-i En 

Rank 1 < 7 ! qf qf -■ 2. (3.13) 

1 <ii qj qj . 


li 



(3.14a) 


Relation (3.13) is, in turn, equivalent to the two scalar equalities 

qiq 2 EN-A - (q \ + q 2 )E N ^ 2 + En-\ - 0, 

q\q 2 Ejss -2 - (<7i + q 2 )Es-\ + E^ = 0, (3.14b) 

which constitute the one-way-to-the-right discrete Helmholtz equation . Relations (3.14a) and (3.14b) 
supplement the scheme (3.2) at n = N - 1 and n = N, respectively, i.e., at the two near-edge nodes of 
the grid where the regular five-point wide stencil of (3.2) cannot be applied. 

The two-way ABC at z = 0 also has to possess the capability of radiation boundary conditions, i.e., it has 
guarantee the transparency of the interface for all the waves that propagate/decay to the left. In other words, 
we require that at the left boundary the outgoing, i.e., scattered, waves be given by E^ c&i) = c\q^ 11 + c 2 q^ n . 
Assuming for a second the homogeneity: g n — 0 near n = 0, we could obtain similarly to (3.13): 

^(scat) ^i(scat) ^i(scat) ^(scat) 

Rank 1 r/j” 1 q^ 2 q f 3 =2. (3.15) 

1 f/.J 1 q^ 2 qC* _ 

Relation (3.15), again, is equivalent to the one-way-to-the-left discrete Helmholtz equation: 

E ( o CM) - (9i + 9-->)d S< ' at> + 9i92^ Sfa ' ’ = 0, (3. 16a) 

E[ scat) - (9i + 9-;)£f at> + 9i9 2 £< scat> = 0. (3.16b) 

Equations (3.16a), (3.16b), however, cannot be immediately used as the ABC at z = 0 because the foregoing 
assumption of homogeneity near the interface is, generally speaking, not correct, and moreover, equations 
(3.16a), (3.16b) do not account for the incoming wave at z = 0 (see Section 2.1), i.e., do not have the 
important two-way capability. The analysis of [12] shows that to accurately address both issues, i.e., the 
inhomogeneity that comes from the previous iteration and the presence of the incoming wave, it is sufficient 
to introduce particular modifications to the right-hand side g n only at two nodes: n — 0 and n = 1. The 
corresponding modification due to the incoming signal is obtained by simply substituting the right-traveling 
incoming wave E^ ttu ^qf into the one-way-to-the-left Helmholtz equation (3.16a), (3.16b). Altogether, the 
two-way ABCs at z = 0 are given by (cf. formulae (3.16a). (3.16b)): 

Eo — (<71 + q 2 )E\ -f q\q 2 E '2 — (Jo, (3.17a) 

E\ ~(q\ + qz)E 2 + qxqiEs = g[, (3.17b) 

where prime denotes the aforementioned modification of the right-hand side, see [12]. Again, relations (3.17a) 
and (3.17b) supplement the scheme (3.2) at the near-edge nodes n = 0 and n = 1, respectively, where the 
regular five-point stencil cannot be applied. Straightforward considerations based on the linear superposition 
principle and uniqueness (see [12]) guarantee that inhomogeneous relations (3.17a), (3.17b) correctly specify 
the incoming signal at z = 0 and still ensure the reflectionless propagation of all the outgoing waves through 
z — 0 toward z — —o o. 

3.4. Computational Complexity. The computational complexity of one solution of equation (3.1) 
is 0(N z N r In N r ) operations, where N z and N r are the corresponding grid dimensions. Indeed, the cost of 
solving each of the Ay one-dimensional systems (3.2) is linear with respect to N z , because each of this systems 
needs to be solved repeatedly for multiple right-hand sides. As such, the sparse LU decomposition can be 
performed only once ahead of time, and the cost of backward substitution is linear. Therefore, the overall 
complexity is dominated by the cost of N z direct and inverse FFTs of length AT r , which is 0{N z N r \nN r ). 


12 


4. Results. In this section we present simulation results for the Gaussian initial conditions Ejj lc — 
exp(— r 2 ) and V’o = («ofcg) 1/4 exp(-r 2 /rg) for the NLH and NLS, respectively. Denoting, as before, the 
input power of the incoming wave by A r (0), we define the fractional input power as 

P = A r (0)/Ay, (4.1) 


i.e., p — 1 when the input power is equal to the NLS critical power A T C « For the Gaussian initial conditions 
used in our simulations p = fco \j2c /Sit [12]. In all simulations we set ko = 8 and ro = 1. 


Table 1.1 

Threshold values of linear damping 


Case No. 

f 

p = a m/s T c 

"t.h 


1 

0.06 

90% 

0 

0 

2 

0.07 

97.5% 

0 

0 

3 

0.072165819 

99% 

0 

0 

4 

3tt/128 

100%. 

9.6 • 10 s 

0 

5 

0.075 

100.9% 

0.00023 

0 

C 

0.08 

104%. 

0.0007 1 

0.00025 

7 

0.1 

116% 

0.0027 

0.0025 

8 

0.125 

s 

130%, 

0.004!) 

0.0062 

9 

0.15 

142% 

0.0070 

0.010 

10 

0.2 

164%, 

0.0145 

0.019 

11 

; 

0.3 

202%, 

0.030 

0.035 

12 

j 

0.4 

233% 

0.044 

0.050 

13 

0.5 

261% 

0.058 

0.065 


In Table 4.1 we show the calculated threshold values ^ and $ t s h . The quantity * ri Table 4.1 represents 
the smallest non-negative value of for which we obtain a global solution of the NLH. By this we mean 
that the nonlinear iterations converge in the sense that the value of max ;tr (£^ ,, + 1 * — E {n) )/ max :(r £ (,,+1) 
drops by at least a factor of 10~ 6 in the course of iterations on the computational domain 0 < z < 40 and 
0 < r < 40, with grid sizes h z = A/20 and h r = A/8, where A - 2i r/fe 0 . The particular choice of the domain 
size and grid resolution is “inherited'’ from our previous numerical experiments, see [10.12]. The values of <5*}, 
in Table 4.1 are obtained with at least two significant digits by repeatedly running the code for a given e and 
varying 5 , which allows one to “close im T on the threshold. As, however, discussed in Section 1, with a larger 
computational domain and/or a finer grid it may be possible* to obtain regular solutions for smaller values 
of <5, hence, to obtain a lower value of the threshold For example, using the same computational domain 
and twice as fine grid: h z = A/40 and h r = A/16, we could obtain <5^ = 0.0133 instead of — 0.0145 
for the data on row No. 10 of Table 4.1 (e = 0.2). Likewise, using the original grid resolution h z = A/20 
and h r = A/8 and the computational domain which was twice as large: - max = 80 and r max = 80. we could 
obtain 6^ = 0.0022 instead of <5^ = 0.0027 for the data on row No. 7 of Table 4.1 (e = 0.1). In other words, 
the values of tfjf, from Table 4.1 should be considered upper hounds for the actual thresholds. However, the 
quantitative limits of pursuing this venue are still unexplored, i.e., it is not known how far down in ^ one 
can go by increasing the domain size and/or grid resolution. Our ability to answer this question is obviously 
limited by computer resources, and as of yet the question remains open . In particular, it is unclear whether 
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we can achieve = 0 by choosing a sufficiently large domain and/or fine grid. 

Similarly, the quantity £ t s h in Table 4.1 represents the smallest non-negative value of damping S for which 
the NLS solution does not blow up. In our NLS simulations we use standard fourth-order finite difference 
schemes for the spatial derivatives and explicit fourth-order Runge-Kutta for marching in z. As has recently 
been shown in [9], in finite-difference simulations of NLS solutions that are known analytically to become 
singular, the computed solution still remains bounded. Therefore, there is always an element of arbitrariness 
in selecting a numerical criterion for blowup in NLS simulations. In our NLH simulations the largest relative 
increase in amplitude due to self- focusing has never exceeded a factor of two. In order to make the blowup 
criteria in NLH arid NLS simulations as close to one another as possible, we define the computed NLS 
solution as becoming singular once its amplitude increases by a factor of two. We checked that altering this 
NLS blowup criterion leads to only minor changes in the results for 5 t s h . For example, using the blowup 
criterion of relative focusing by a factor of 4, rather than 2, for e = 0.08 (row No. 6 of Table 4.1) gives 
<5 t s h = 0.00021 instead of <5 t s h = 0.00025; and using this new criterion for e = 0.15 (row No. 9 of Table 4.1) 
yields — 0.0089 instead of 6f h = 0.010. In particular, this change does not affect our main finding of 
initial conditions for which <5^ < Sf h . 

As expected, for both the NLS and the NLH the threshold values of 6 increase with e (i.e., a larger 
amount of damping is needed to arrest collapse of beams with higher input power). For e = 0.06 and 
e = 0.07 the input power is below critical. Therefore, both the NLS and the NLH have global solutions 
for (5 = 0. This behavior for the NLH holds (at least.) till e = 0.072165819, which corresponds to the last 
subcritical value 5 that we have checked: A r (0) being equal to 99% of N c . 

Starting from e = 37r/2fco « 0.073631077, which corresponds exactly to A r (0) = N c , the NLH requires a 
certain positive amount of damping S to maintain the regularity of the solution. For the NLS, the solution 
with no damping remains regular till e — 0.75, which corresponds to p — A r (0)/A T c = 1.009. Indeed, it is 
known that N v is only a. lower bound for the threshold power for NLS collapse, and that any initial condition 
which does blow up, and whose amplitude |^q| is not equal to the ground state profile 3 1 ^ 4 y^sech(2r), has 
power strictly above N v [7, 16, 17]. In our simulations we have discovered that for e = 37t/2/lq, which is the 
critical value for the NLH, as well as for the moderately supercritical values e — 0.075, e = 0.08 and e = 0.1, 
when the input power A r (0) is only slightly above N v , the threshold damping for the NLH is larger 6 than that 
for the NLS: (5 t ^ > Sf h . However, for input powers that are equal or higher than 1.30 A/, (which corresponds 
to € = 0.125) this trend reverses, see Table 4.1, and we obtain <5^ < Sf h . Thus, for A%0) > 1.30A r c , 7 there 
must be other mechanisms in the NLH not present in the NLS that help suppress the formation of singularity 
in the solution. Therefore, we may conclude that in this regime nonpamxiality and backscattering help arrest 
collapse of nonlinear waves. 

In [6] Fibich has used asymptotic analysis to show that 

&~c(p-l) 8/2 , (4.2) 

where p is the fractional critical power (4.1). In Figure 4.1 we put this theoretical prediction to a test by 
plotting the values of and <5^ as a function of (p — 1). When we computed the best fit of the values of 
<5^ h with the two-parameter family of curves <5 t ^ = c(p — 1) Q , we obtained a = 1.517, which is in excellent 

5 As mentioned in Section 1 . for larger subcritical values of A T ( 0) the convergence of nonlinear iterations becomes prohibitively 
slow. 

6 We remind that the values of 5*^ in Table 4.1 are only upper bounds for the threshold; lower values may be obtained by 
refining the grid and/or enlarging the computational domain. 

7 More precisely, A r (0) higher than some value between 1.16A' C and 1.30A ; C 
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agreement with formula (4.2). Relation (4.2) also provides a good approximation to the data points <5^, 
see Figure 4.1. The only exception is the lowest- power NLH data point in Figure 4.1 that corresponds to 
f = 0.08 (row No. G in Table 4.1), for which the value of S\l has most likely been overprodicted numerically 
because of the computational constraints discussed previously. 



Fig. 4.1. Threshold values ( hollow bullets “o ”) and (asterisks “*”) as a f motion of {p- l) for the data in Table 4-T 
The solid line 0.035(p — l) 1 517 is the best Jit to the values of Sf h . 

In Figure 4.2 we plot the on-axis (r = 0) amplitudes of the NLH and NLS solutions for e = 0.2 and 
various values of <5. The on-axis behavior is most representative of the ph\ sical processes that we are studying, 
because for symmetric beams this is the location of the peak intensity. When ^ = 0.0145, the NLH 

solution exists globally but the NLS solution becomes singular at a finite propagation distance. As the value 
of damping increases, both the NLS and the NLH solutions undergo less focusing. For all the cases for 
which both solutions remain regular, the NLS solution curve is higher i han the NLH one from z = 0 until 
its maximum, i.e., the point of the arrest of collapse. This provides additional support to the foregoing 
conclusion that nonparaxiality and backscattering arrest collapse of nonlinear waves. Note that after the 
collapse has been arrested, the NLS solution becomes lower than the NLS one. One possible explanation for 
this is that the NLS solution is undergoing higher focusing, hence it losos more power due to damping. 

We emphasize that at z — 0 the NLH solution is not equal to £• J 1C , see Figure 4.2. The difference 
between the two is due to backscattering, and can be used to quantify the level of backscattering for a 
particular setting, see [10, 12]. 8 In Table 4.2 we provide the values of maximum self-focusing and maximum 
backscattering in the NLH, defined as max r , 2 \E(z,r)\ and max r |£(0, r) — £f] ir (r)|, respectively, for various 
values of e and S . The dash “ in a particular cell of Table 4.2 means that the level of damping was 

H Thcn* are, in fact., two phenomena that account for the discrepancy between the NLH and NLS curves: Nonparaxiality of 
t he forward propagating wave and backscattering. Because the problem is nonlinear, these two mechanisms cannot be easily 
and explicitly told apart inside the domain. The only location where we can rlearlv say that the difference is purely due to 
backscattering is the “inflow” interface 2 = 0, see [10]. 
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8 = 0.0145 



6 = 0.019 



8 = 0.025 


8 = 0.040 




FIG. 4.2. On-axis amplitude of NLH (solid) and NLS (dashes) solutions for e = 0.2 and various values of 6. 


insufficient to guarantee the convergence of the numerical algorithm. As expected, for a given level of 
damping £, the NLH solution undergoes stronger self-focusing as the nonlinearity coefficient e increases. The 
level of hackscattering also increases with the increase of e . As also expected, for a given input power e, 
when the damping S increases the NLH solution undergoes weaker self- focusing (see Figure 4.2). Surprisingly, 
however, changing the value of damping S has very little or no effect on the level of backs cattering. To further 
corroborate this observation, we picked a particular value of the nonlinearity coefficient, e = 0.2, and ran an 
additional series of numerical tests with a substantially more fine sampling for 6. These results, which are 
presented in Table 4.3, confirm that backscattering is not affected by linear damping. This phenomenon 
certainly cannot be explained by saying that linear damping has the overall negligible effect, since its effect 


Table 4.2 

Maximum absolute levels of self-focusing and backscattering in the NLH for a variety of e and S. 



Maximum self-focusing 

Maximum backscattering 

6 = 0.0145 

S = 0.0175 

<5 = 0.0210 

S = 0.0145 

S = 0.0175 

5 = 0.0210 

e = 0.15 

1.1179 

1.0601 

1.0162 

0.0372 

0.0373 

0.0373 

— 1K1 






0.0421 

B319i 



1.1716 


0.0466 


e = 0.225 

— 

— 

1.3242 

— 

— 

0.0509 
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on the focusing dynamics can be clearly seen through the changing values of the maximum focusing both in 
Table 4.2 and in Figure 4.2. At present., we have no good explanation to this surprising observation. 


Table 1.3 

Maximum absolute levels of self-focusing and backs cattcring in 'he NIAl for f = 0.2. 


Case No. 

Damping S 

Max. self-focusing 

Max. backscattering 

1 

0.0145 

1.5515 

0.0465 

2 

0.0147 

1.5296 

0.0465 

3 

0.0150 

1.4992 

0.0465 

4 

0.0155 

1.4538 

0.0465 

5 

0.0160 

1.4135 

0.0466 

6 

0.0165 

1.3776 

0.0466 

7 

0.0170 

1.3451 

0.0466 

8 

0.0175 

1.3158 

0.0466 

9 

0.0180 

1.2892 

0.0466 

10 

0.0190 

1.2428 

0.0466 

11 

0.0200 

1.2041 

0.0466 

12 

0.0210 

1.1716 

0.0466 


5. Concluding Remarks. The question whether nonparaxiality and backseat tering may arrest 
collapse of nonlinear waves has been open for many years. While the answer to this question is probably 
positive, no conclusive argument toward it, whether analytical or numerical, has been previously available in 
the literature. In this study we addressed this question within the framework of the linearly damped NLH 
and NLS. As has been mentioned, the addition of linear damping is not ad-hoc because it has both physical 
and mathematical motivation. Methodologically, linear damping provides a very useful ’‘extra dimension 
that allowed us to efficiently control the solutions of the NLH arid NLS. Specifically, the variation along this 
extra dimension has helped us to numerically identify the regimes, for which the NLS solution blows up, 
while the NLH solution remains regular. In other words, our results furnish the first ever definite numerical 
evidence that nonparaxiality and backscattering can arrest collapse . The question whether regular solutions 
to the NLH still exist in the absence of damping remains open as of yet. However, we hope that the arguments 
based on linear damping and the limiting absorption principle may be useful for proving global existence 
and uniqueness, both for the damped NLH and for the undamped NLH 
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